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Abstract — In orthogonal frequency-division multiplexing (OFDM) 
systems operating over rapidly time-varying channels, the orthogonality 
between subcarriers is destroyed leading to inter-carrier interference 
(ICI) and resulting in an irreducible error floor. In this paper, a new 
and low-complexity maximum a posteriori probability (MAP) detection 
algorithm is proposed for OFDM systems operating over rapidly 
time-varying multipath channels. The detection algorithm exploits 
the banded structure of the frequency-domain channel matrix whose 
bandwidth is a parameter to be adjusted according to the speed of the 
mobile terminal. Based on this assumption, the received signal vector 
is decomposed into reduced dimensional sub-observations in such a 
way that all components of the observation vector contributing to the 
symbol to be detected are included in the decomposed observation 
model. The data symbols are then detected by the MAP algorithm 
by means of a Markov chain Monte Carlo (MCMC) technique in an 
optimal and computationally efficient way. Computational complexity 
investigation as well as simulation results indicate that this algorithm 
has significant performance and complexity advantages over existing 
suboptimal detection and equalization algorithms proposed earlier in 
the literature. 

Index Terms- OFDM, MAP detection, Monte carlo technique, Gibbs 
sampling. Intercarrier interference, fast time-varying channels. 



I. Introduction 

Orthogonal frequency-division multiplexing (OFDM) has been 
shown to be an effective method to overcome inter-symbol inter- 
ference (ISI) caused by frequency-selective fading with a simple 
transceiver structure, and is consequently used in several existing 
wireless local and metropolitan area standards such as the IEEE 
802.11 and IEEE 802.16 families. IEEE 802.11 wireless LAN 
(WLAN) technology has become very popular for providing data 
services to Internet users although its overall design and feature set 
are not well suited for outdoor broadband wireless access (BWA) 
applications [1]. Therefore, IEEE 802.16 has been developed as a new 
standard for BWA applications [2]. Recently, the much-anticipated 
Worldwide Interoperability for Microwave Access (WiMAX) tech- 
nology was introduced to promote the 802.16 standards while in- 
troducing features to enable mobile broadband services at vehicular 
speeds beyond 120 km/h. 

OFDM eliminates ISI and simply uses a one-tap equalizer to 
compensate multiplicative channel distortion in quasi-static channels. 
However, in fading channels with very high mobility, the time 
variation of the channel over an OFDM symbol period results 
in a loss of subchannel orthogonality which leads to inter-carrier 
interference (ICI). A considerable amount of research on OFDM 
receivers for quasi-static fading has been conducted, but a major 
hindrance to such receivers is the lack of mobility support [3], Since 
mobility support is widely considered to be one of the key features 

E. Panayirci is on sabbatical leave from Kadir Has University, Istanbul, Turkey. 
This research has been conducted within the NEWCOM++ Network of 
Excellence in Wireless Communications and WIMAGIC Strep projects funded 
through the EC 7th Framework Programs and was supported in part by the 
U.S. National Science Foundation under Grant CNS-0625637. 



in wireless communication systems, and in this case ICI degrades 
the performance of OFDM systems, OFDM transmission over very 
rapidly time varying multipath fading channels has been considered 
recently in a number of recent works [4-15]. 

The techniques proposed in these works range from linear equal- 
izers, based on the zero-forcing (ZF) or the minimum mean-squared 
error (MMSE) criterion [4—13], to nonlinear equalizers based on 
decision-feedback or ICI cancelation [9-14]. Also near maximum- 
likelihood approaches have been proposed [16]. It has been shown 
that nonlinear equalizers based on ICI cancelation generally out- 
perform linear approaches [10-13]. However, linear equalizers still 
preserve their importance mainly because they are less complex. 

In [12], performances of Matched Filter (MF), Least Squares (LS), 
Minimum Mean Square Error (MMSE) and MMSE with Successive 
Detection (SD) techniques with optimal ordering have been investi- 
gated. However, since the number of subcarriers is usually very large 
in high speed wide-band wireless standards, even the linear MMSE 
equalizer considered in [12] demands very high computational load. 

The specific structure of the Doppler-induced ICI in OFDM 
systems operating over highly mobile channels presents a distinctive 
feature of limited support of the Doppler spread that can be exploited 
by the receiver. References [6-10] exploit the banded character of the 
frequency-domain channel matrix to reach a complexity that is only 
linear in the number of subcarriers. In a certain sense, the assumption 
of a banded frequency-domain channel matrix is a natural extension 
of the time-invariant channel case, in which the frequency-domain 
channel matrix is diagonal and hence banded with the smallest 
possible bandwidth. 

In [7], using the banded structure of the channel matrix, a simple 
frequency domain equalizer has been proposed that can compensate 
for the loss of subchannel orthogonality due to ICI. However, the 
detection performance of the technique degrades substantially, since 
the data to be detected cannot fully use the contributing observation 
elements. The work presented in [13] combined [7] and [10] to derive 
a recursive decision feedback equalizer receiver for ICI suppression. 
The iterative MMSE serial linear equalizer (SLE) of [10], which takes 
the banded structure of the channel matrix into account, seems to be 
one of the most promising approaches to compensate for ICI. Iterative 
MMSE is then applied to estimate frequency-domain symbols. In [8], 
a block MMSE equalizer for OFDM systems operating over time- 
varying channels is presented. By exploiting the banded structure of 
the frequency-domain channel matrix, the complexity of the resulting 
algorithm turns out to be smaller than that of [10]. 

In this paper, a new computationally feasible, maximum a pos- 
teriori probability (MAP)-based data symbol detection algorithm is 
proposed for OFDM systems operating over highly mobile channels, 
as an alternative to the existing suboptimal equalization/detection 
techniques summarized in the above paragraphs. The proposed de- 
tection algorithm exploits the banded structure of the frequency- 
domain channel matrix whose bandwidth is a parameter to be adjusted 
according to the speed of the mobile terminal. This assumption 
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enables us to decompose the main received signal vector into finite 
numbers of reduced-dimensional, sub-received signal vectors from 
which the data symbols can be detected by the MAP algorithm in 
an optimal and computationally efficient way. The decomposition 
is achieved in such a way that all the components of the received 
vector that contribute to the symbol to be detected are included in the 
decomposed observation model. Data symbols in each sub-received 
signal model are then detected successively by a MAP detection 
algorithm. To implement MAP symbol detection in a computationally 
efficient way, we employ a Markov Chain Monte Carlo (MCMC) 
technique based on Gibbs sampling, which is a powerful statistical 
signal processing tool to estimate a posteriori probability (APP) 
values. 

The resulting detection algorithm is compared with previously 
proposed algorithms in terms of both bit error rate (BER) and 
complexity requirements. Computational complexity investigation as 
well as simulation results indicate that our algorithm has significant 
performance and complexity advantages over the existing suboptimal 
detection and equalization algorithms. 

II. System Model 

Let us consider an OFDM system with N subcarriers and available 
bandwidth B — 1/T„ where T 3 is the sampling period. A given 
sampling period is divided into N subchannels by equal frequency 
spacing A/ = B/N. At the transmitter, information symbols are 
mapped into possibly complex-valued transmitted symbols according 
to the modulation format employed. The symbols are processed 
by an N— length Inverse Fast Fourier Transform (IFFT) block that 
transforms the data symbol sequence into the time domain. The time- 
domain signal is extended by a guard interval containing G samples 
whose length is chosen to be longer than the expected delay spread 
to avoid ISI. The guard interval includes a cyclically extended part 
of the OFDM block to avoid ICI. Hence, the complete OFDM block 
duration is P — N + G samples. The resulting signal is converted to 
an analog signal by a digital-to-analog (D/A) converter. After shaping 
with a low-pass filter (e.g. a raised-cosine filter) with bandwidth B, 
it is transmitted through the transmit antenna with the overall symbol 
duration of PT S . 

Let h(m, I) represent the Ith path (multipath component) of the 
time-varying channel impulse response at time instant t = mT s . The 
discrete-time received signal can then be expressed as follows: 

L-l 

y(m) = h(m, l)d(m — I) + w(m), (1) 
;=o 

where the transmitted signal d(m) at discrete sampling time m is 
given by 

1 N ~ 1 

d(m) = ^=J2 s ( k y 2 ™ k/N > <?) 

L is the total number of paths of the frequency selective fading 
channel, and w(m) is additive white Gaussian noise (AWGN) 
with zero mean and variance _E{|u>(m)| 2 } = o~^. The sequence 
s(k),k = 0, 1, ■ • • , N — 1, in (2) represent either quadrature- 
amplitude modulation (QAM) or phase-sift-keying (PSK) modulated 
data symbols with _E{|s(fc)| 2 } = 1. 

At the receiver, after passing through the analog-to-digital (A/D) 
converter and removing the cyclic prefix (CP), a fast Fourier trans- 
form (FFT) is used to transform the data back into the frequency 
domain. Lastly, the binary data is obtained after demodulation and 
channel decoding. 

The fading channel coefficients h(m, I) can be modeled as zero- 
mean complex Gaussian random variables. Based on the wide- 
sense stationary uncorrelated scattering (WSSUS) assumption, the 



fading channel coefficients in different paths are uncorrelated with 
each other. However, these coefficients are correlated within each 
individual path and have a Jakes Doppler power spectral density 
having an autocorrelation function given by 



E{h{m, l)h*(n, I)} = a 2 hl J (2nf d T s {m - n)), 



(3) 



where af Z[ denotes the power of the channel coefficients of the 
ith path. f d is the Doppler frequency in Hertz so that the term 
fdT s represents the normalized Doppler frequency of the channel 
coefficients. Jo(.) is the zeroth order Bessel function of the first kind. 
By using 10 in 10, the received signal can be written as 

JV-l L-l a fc( _ () 

y(m) = -= J2 s{k)^h{m,iy^ k « ' +w(m), (4) 



(5) 



fc=0 1=0 

which upon defining the time-varying channel transfer function 



H(k,m) = ^ h(m,l)e- j2nl 



fc/JV 



becomes 



y(m) = -L V s(k)H(k, m)e j27rmfc/JV + w(m). (6) 

The FFT output at the k th subcarrier, after excluding the guard 
interval, can be expressed as 

» m= 

= s(k)G(k,k)+I(k) + W{k), (7) 

where I(k) is ICI caused by the time-varying nature of the channel 
given as 



N-l 

i(k)= J2 «(0<w) 



(8) 



G(k, i) in ^ represents the average frequency domain time- 
varying channel response, defined as 



G(k,i) = (1/N) H(i,m)e 



j2irm(i-k)/N 



(9) 



Similarly, the term G(k, k) = i 5Zm=o H{k,m) in l|9} represent 
the portion of the average frequency domain channel response at the 
kth subcarrier and W(k) denotes discrete Fourier transform of the 
white Gaussian noise w(m): 



-j27vmk/N 



(10) 



Because of the term I(k) in 0, there is an irreducible error floor 
even in the training sequences since pilot symbols are also corrupted 
by ICI, arising from the fact that the time-varying channel destroys 
the orthogonality between subcarriers. Therefore, channel estimation 
should be performed either jointly with data or before the FFT block 
in order to compensate for the ICI. Note that if the channel is 
static or very slowly time-varying, that is H(k,m) ~ H(k), then 
it can be easily shown that G(k, k) = H(k) and I(k) — for 
k = 0, 1, • • • ,N — 1, resulting in a received signal at the output of 
the FFT processor corresponding to the fcth OFDM symbol given by 

Y{k) =H(k)s(k) + W(k). 

From (6) and (7), the FFT output received signal can be expressed 
in vector form as 

Y = Gs + W (11) 
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where Y = [Y(0), Y(l), Y(N- 1)] T , s = [s(0), s(l), .., s(iV- 
1)] T and W = [W(0), W(l),..., W(N - 1)] T . For k,i = 
0, 1, • • • , JV - 1, the (k, i)th element of the matrix G = [G(k, i)] G 
qNxn re p resen ti n g the time-varying channel is given by (9). 

Under the assumption that the channel matrix G in (11) is per- 
fectly known at the receiver, the maximum likelihood (ML) detector 
performs an exhaustive search over the entire set of signal vectors 
whose components are selected from the signal constellation formed 
by the modulation scheme chosen. Especially in IEEE 802.16 based 
systems, the length TV of each OFDM symbol is very large; it can 
take values as large as N = 1024 or even N — 2048 especially for 
high mobility applications. In this case an exhaustive search for the 
ML solution would be very complex since the search space has an 
extremely large number of constellation points, (jSI^, where |S| is 
the cardinality of the signal constellation). On the other hand, all of 
the lower complexity linear detectors given in Table 1 are suboptimal 
since they do not take into account the correlation of the components 
of the transformed noise and yield noise enhancement. Recently, a 
nonlinear recursive detection technique using the decision-feedback 
principle, namely the MMSE-SD algorithm (VBLAST), has been 
proposed [12], The performance of VBLAST depends critically on 
the order in which the data vector components are processed. To 
minimize error propagation effects and to improve the detection of 
unreliable components, more reliable data vector components should 
be detected first. Therefore, the algorithm depends on calculation 
of the post-detection signal-to-interference-plus-noise ratio (SINR) 
based upon MMSE detection as a measure of reliability, and so the 
calculation of SINR is compulsory at each iteration. Therefore this 
algorithm is computationally intensive as a number of pseudo inverse 
operations need to be performed. Moreover, its complexity grows 
exponentially with the total number of subcarriers. 



TABLE I 
Linear Detection Methods 



Method 


Solution 


Matched Filter (MF) 


g = Q{G T Y} 


Zero Forcing (ZF) 


§ = q{(gtg)- 1 gty} 


MMSE 


s = Q{(GTG + I N ^)- 1 GTY} 



where Ijv is the N-by-N identity matrix. 

III. MAP Detection Algorithm 

The necessity of detecting large numbers of symbols in OFDM 
systems employed especially in highly mobile and wide-band wireless 
systems represents a significant computational burden as well as 
creating some convergence problems. However, as is known [9], 
time-varying channels produce a nearly-banded channel matrix G 
whose only main diagonal, Q subdiagonals and Q superdiagonals are 
nonzero. The bandwidth 2Q + 1, which is defined as the total number 
of non-zero diagonals in G, is a parameter to be adjusted according 
to the mobility-rate of the channel. The significant coefficients of 
G thus are confined to the 2Q + 1 central diagonals. The parameter 
Q £ {0, 1, • • • N/2— 1} controls the target ICI response length; larger 
Q corresponds to a longer ICI span and, thus, increased estimation 
complexity. In general, Q should be chosen proportional to the width 
of the Doppler spectrum of the channel. 

We now present an optimal low-complexity MAP detection algo- 
rithm to detect the data symbols s from Y taking into account the 
banded structure of the channel matrix G. From the observation Y 
in (11), the receiver attempt to detect the OFDM output symbol s, 
assuming that G is completely known by the receiver. The banded 
structure of the channel matrix G implies that the data symbol 



s(k), k = 0, 1, ■ • ■ , N — 1, contributes to a maximum of 2Q + 1 
observation elements as follows: 

Y fc = [Y(j k ),Y(jk + l),--- ,Y(i h )f, for (12) 

jk = max{0, k — Q} and i k = min{iV — 1, k + Q}. 

Based on this observation, the received signal in (11) can be de- 
composed into N reduced-dimensional sub-observations from which 
the data symbols can be detected in an optimal and computationally 
efficient way. For a given index k = 0, 1, • • ■ , N — 1 and Q, it can 
be easily shown from (11) and (12) that 

Y fc = G fcSfc + W fc (13) 
where, s fe = [s(j k )), W k = [W(j k )], G k = [G(i k ,j k )], for 

I Lk ± max{0, k-Q} <i k < I Uk ± min{7V - 1, k + Q} (14) 
and 

J Lk = max{0,fc-2Q} < j k < J Uk = min{iV- 1, k + 2Q}. (15) 

Note that due to the banded structure of G, some elements of the 
matrices G fe are zero and dim(G fe ) < (2Q + 1) x (2(2Q + 1) - 1). 
The Gfe's reach their maximum dimension when 2Q + 1 < k < 
(N-(2Q + 1)). 

For k = 0, 1, • • • , N — 1, the MAP estimate of the data symbol 
s(k) given Yfc is 

'sMAp(k) = s(k) = arg max P(s(k)\Y k ), (16) 

s(fc)£S 

where 5 denotes the set of signal constellation points from which 
s(k) takes values. Based on this approach, s(k) can be detected 
sequentially for k = 0, 1, • • ■ , N — 1, incorporating the outcomes of 
the previous estimates in a decision-feedback mode as follows. 

■ For k = 0, determine the estimate S"(0) from (16). 

■ For k = k + 1 modify the observation vector Yfc by subtracting 
the terms coming from the contributions of the estimated data 
symbols s"(0), s"(l), • • • ,~s(k — 1) as 

fc-i 

Yfc 4 Yfc - J2 S k l) s(l) (17) 

where g?' is the Zth column of Gfc and JL k and Jjj k are defined in 
(15). 

■ Determine the MAP estimate of *s(k) from 

Yfc = Gfci fc + Wfc (18) 

as 

s(k) = arg max P(s(k)\Y k ), (19) 

where 

G^Gfc-[ g r,gW,..., g r 1 )], (20) 
and Sfc is the vector obtained by removing the first k — 1 elements 

Of Sfc. 

■ END IF k = N - 1. 

The major problem is finding the values of SMAp(k) in a com- 
putationally efficient manner. To see this difficulty we assume that 
the data symbols are independent and identically distributed binary 
phase shift keying (BPSK), taking values of +1 and —1. Note that 
higher dimensional signal constellations can be treated similarly with 
a straightforward extension. The conditional probability of H(k) given 
the observation vector Yk can be expressed as 
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P(s(k) = +l\Y k ) = J2P(s(k) = +l,%|Y fc ) 
5 r 

= J^P(s(k) = +l|%,Y fe )P(%|Y fc ) (21) 
s r 

where the second identity follows by applying the chain rule of 
probability. The vector in (21) is obtained by canceling the 
component s(k) in Sk and the summation is over all possible values 
of s^-. The number of combinations that s^ takes grows exponentially 
with the dimension of Sk and thus becomes prohibitive for large 
values of the size of this vector. Thus, we resort to the Gibbs sampler, 
a Monte Carlo method to calculate the a posteriori probabilities of 
the unknown symbols. 

A. MAP Detection Based on Gibbs Sampling 

The Gibbs sampler is an MCMC sampling method for numerical 
evaluation of multidimensional integrals. Its popularity is gained from 
the facts that it is capable of carrying out many complex Bayesian 
computations. In this section we briefly explain the application of the 
Gibbs sampling technique to our symbol detection problem where the 
observation process is given by (18). For notational convenience we 
drop the index k and the "tilde" from all the involved variables, e.g. 
Y, G are shorthand notations for Yfc, Gfc, respectively. (21) can then 
be expressed as 

P(a(*) = +1|Y) = ]rP( S (fc) = +l|s ¥) Y)P(s ¥ |Y) 

= £ S¥ | Y {P( S (fc) = +l|Sfc,Y)}. (22) 

According to the Gibbs sampling based statistical Monte Carlo 
estimation technique, an estimate of (22) can be evaluated by taking 
the empirical average 



1 " 

P(s(k) = +1|Y) = JL P(s(k) = +l|si"', Y) 



(23) 



where si. for n = 1,2, ••• ,N S are samples drawn from the 
conditional distribution P(sj|Y), There is a substantial body of 
literature concerning the Monte Carlo Gibbs sampling technique; see, 
e.g., [17], [18]. One possible version of the Gibbs sampler suitable for 
calculating the a posteriori probabilities in (21) may be summarized 
as follows. 

Let s = [s(0), s(l), • • ■ , s(N — 1)] T be a vector of unknown data 
symbols. Let Y be the observed signal. To generate random samples 
from the distribution P(s|Y), given the samples from the (n — l)th 
iterations^" 1 ) = [s^ 1 ) (0), s^" 1 ' (1), • • • , (N - 1)] T , the 

Gibbs algorithm iterates at the nth iteration as follows to generate 



(N-l)] 
,N-1, 

(n) (n-l) 



n) \ 



the samples s (n) = [s (n) (0), s (n) (1), 

■ Initialize s^ ' randomly; 

■ for n = 1, 2, • • • ,N T and for k = 0, 1 • 
draw sample from P (s(k)\sg^ 

Note that to ensure convergence, the Gibbs iteration is usually 
carried out for Nt = Nb + N s iterations. The first Nt iterations 
of the loop is called the burn-in period which is necessary for the 
Monte Carlo simulation to reach its stationary distribution. Only the 
samples s (n) = [s { n) , s[ n) , ■ ■ ■ , 4?li] T , n = N b + 1, • ■ • , N T , from 
the last N s iterations are used to calculate the expectation. 

It is known that under regularity conditions [18, 19], 

(i) the distribution of s' n ' converges to P(s|y), as n — > oo. 

(ii) (V^E^iW) = +l|4 n) ,Y) = J2 s _P(s(k) = 
+l|Sfc, Y)P(s r |Y), as n -> oo. 



B. Implementation of the Symbol Detector 

From the previous section, we recall that to compute P(s(k) = 
+ l|Y fc ) = Y2 s - P ( s ( k ) = +Mfcl Y fc)> we need t0 perform the 
summation on the right-hand side of (22). When s has a large 
dimension, the exact evaluation of this summation may not be 
feasible and other more efficient techniques must be adopted. In this 
section the Gibbs sampling-based Monte Carlo method summarized 
in the previous section will be applied to develop a computationally 
efficient algorithm for calculation of the a posteriori probabilities 
P(s(fc)]Y). From (23), it follows that we need to evaluate P(s(k) = 



+l|si n) ,Y), forn = 1,2, ■•■ , Nt- For this we define 

, , P(s(k) = +l|si n) ,Y) 
A " ' 4 In V W 



P(s{k) = 
from which it can be easily seen that 

P(s(k) = +l\4 n \Y) = ■ 



l|4 n) ,Y) 



exp 



(24) 



(25) 



A^.™' can be computed by expanding P(s(k) = +l\s^ l> ,Y) as 



(«(*) 



= +l|s^ n) ,Y 



(26) 



p(Y| s (fc) = +l,s^ ) )P( S (fc) = +l,s. 



e. W6{+ i,-i } p (yi-(*) = +i, 4 n) ) p (s(k) = +i, 4" 

The data symbols are assumed to be independent and equally likely, 
Therefore, it follows from (26) and (24) that 



p(Y| S (fc)=+l )S ( n) ) 



(27) 



P (YKfc)=-l,4"> )' 

Since p(Y|s) ~ exp(— |y — Gs| 2 ), after some algebra, (27) can be 
expressed as 



A"( Sfe ) = ^5R{gI(Y-G ¥ s ¥ )} 



(28) 



where denotes the conjugate transpose and 5R{-} denotes the 
real part of its argument. G^ is G with its fcth column gk removed. 
In summary, for k = 0, 1, • • • , N — 1, to estimate the a posteriori 
probabilities P(s(fe)|Y) in (21), the Gibbs sampler runs over all 

( In) ) Nt 

symbols N s times to generate a collection of vectors < sL' > 

I k ) n=N b +l 
which are used in (23) to compute the desired quantities. 

C. Complexity Requirements 

The computational complexity of the MAP symbol detector based 
on Gibbs sampling proposed in this work is determined by the 
parameters N S ,Q,N and the constellation size of the transmitted data 
symbols. The computation of Yfc in (18) for k = 0, 1, • • • , N — 1 
requires a maximum of (4Q 2 +2Q)N complex multiplications (CMs) 
and (4:Q 2 +2Q)N complex additions (CAs) per data block. Assuming 
BPSK signaling, the computation of the a posteriori probabilities 
in (22) requires a maximum of (8Q 2 + 2Q + 1)NN S CMs and 
(8Q 2 + 2Q— 1)NN S CAs and computation of the empirical average 
of a posteriori probabilities of the data symbols in (23) requires AW S 
CSs. Therefore, the whole algorithm requires of maximum N(4Q 2 + 
2Q + (4Q 2 + 2Q)Ag CMs and N(4Q 2 + 2Q + (4Q 2 + 2Q)N S ) 
CAs, leading to a total of 27Y(4Q 2 + 2Q + (4Q 2 + 2Q)N S ) complex 
operations. 

Several low complexity equalization algorithms have been de- 
veloped recently, of which several are worth mentioning here to 
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Fig. 1. BER comparison of various detection algorithms for OFDM systems; 

TU Channel (6 taps), 420 km/h Pig- 2- BER comparison of various detection algorithms for OFDM systems; 

TU Channel (6 taps), 120 km/h 



compare their computational complexities with that of the Gibbs- 
based algorithm. 

Ruguni et al. [8] proposed a block MMSE technique based on 
exploiting the banded structure of the channel matrix G. The matrix 
inversion was obtained using a low-complexity decomposition such 
as Cholesky or the LDL* decomposition. The algorithm requires 
a total of (8Q 2 + 22Q + 4)iV complex operations. Schniter [10] 
proposed a linear serial equalizer also based on exploiting the banded 
structure of the channel matrix. This algorithm requires a total of 
(8/3Q 3 + 2Q 2 + 5/3Q + 4)N complex operations. The complexity 
of the serial MMSE equalizer is higher than that of the block MMSE 
equalizer. 

In the VBLAST algorithm, matrix inversion is needed of dimension 
equal to the number of OFDM subcarriers. As a result, the compu- 
tational complexity of the VBLAST receiver increases rapidly with 
the number of subcarriers, which makes its real-time implementation 
prohibitive for large numbers of subcarriers. 

As can be seen easily, the complexity of our algorithm is of the 
same order of the above equalization algorithms and is lower than the 
VBLAST algorithm. However, as remarked earlier, these algorithms 
are suboptimal as opposed to our optimal MAP detection algorithm 
and perform poorly especially when the ICI is high. It is also worth 
mentioning that our algorithm can be easily extended to an iterative 
multiuser MAP detection scheme for OFDM systems. 



> MMSE 

C> VBLAST 

> MAP GS Q=6. 
MAP Exact Q= 




Fig. 3. BER comparison of various detection algorithms for OFDM systems; 
BU Channel (6 taps), 420 km/h 



- SNR-20dB, V=120 kmh 

- SNR-25dB, V=120 kmh 

- SNR-20dB, V=420 kmh 

- SNR-25dB, V=420 kmh 



IV. Simulation Results 

This section presents computer simulation results of the proposed 
detection methods for rapidly varying mobile radio channels. The 
system operates with a 5 MHz bandwidth and is divided into 512 
tones (N = 512) with a total symbol period (T a ) of 115 /xs, of 
which 12.8 /is constitute the CP. One OFDM symbol thus consists 
of 576 samples, sixty-four of which constitute the CP. The normalized 
Doppler frequencies are fdi *T S — 0.0307 and fda * T s = 0.1075, 
corresponding to an IEEE 802. 16e mobile terminal moving with 
speeds of 120 km/h and 420 km/h, respectively, for a carrier 
frequency of 2.4 GHz. The wireless channels between the mobile 
antenna and the receiver antenna are modeled based on a realistic 
channel model determined by the COST-207 project in which Typical 
Urban (TU) and Bad Urban (BU) channel models are considered. For 
each OFDM symbol, Gibbs sampling is performed for 30 iterations, 
with the first 10 iterations as the burn-in period. 




Fig. 4. BER performance of the proposed MAP algorithm based on Gibbs 
Sampling for the TU channel as a function of Q values 
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The information symbols are BPSK modulated to yield 5 = {±1}. 
Figs. 1-3 compare the BER performance of the proposed Gibbs- 
based MAP detection algorithm, an equalization technique based on 
zero forcing (ZF) proposed in [7], linear MMSE, and the VBLAST 
algorithm proposed in [12], as a function of energy per bit to noise 
power ratio (Eh /No) where JVo is equal to a^. 

ZF causes noise enhancement while it eliminates the ICI. There- 
fore, it is seen that ZF performance is the worst. The reason for this 
can also be explained by the ill-conditioned matrix (G'G) to be 
inverted. This problem can be solved with MMSE, which provides 
a good trade-off between ICI cancellation and noise elevation by 
using the knowledge of the noise level. In other words, the noise 
enhancement can be reduced by the insertion of the noise power in 
the inverse matrix that given in Table 1. 

We note that linear equalization of the received signal is suboptimal 
as mentioned in Section II and hence Gibbs-based MAP detection 
algorithm (MAP-GS) is proposed in Section III. It is observed that 
MAP-GS outperforms both ZF and MMSE receivers while it has 
similar performance to VBLAST. Moreover, the exact MAP perfor- 
mance is also included to benchmark the proposed algorithms. It can 
be concluded from these figures that the compared algorithms have 
similar performance for the speed of 120 km/h while the performance 
difference is obvious for the speed of 420 km/h. Moreover, it is seen 
that the BER performance of all algorithms has slightly decreased for 
the BU channel. In particular, it is observed that a savings in about 2 
dB is obtained at BER = 10 -3 , as compared with MMSE detection 
for the TU channel. 

It has been shown that when a proper detection technique is 
adopted, the time-varying nature of the channel can be exploited 
as a provider of time diversity [12]. In [12], it was demonstrated 
that VBLAST fully utilizes the time diversity while suppressing 
the residual interference and the noise enhancement. Similarly to 
VBLAST, we have seen that MAP-GS is also a useful detection 
technique for time-varying channels while having lower complexity. 
Therefore, in particular, it is not surprising that in simulations the 
performance at 420km/h is better than that at 120km/h. 

Finally, the BER performance of the proposed algorithm is pre- 
sented as a function of Q in Fig. 4. The parameter Q can be chosen 
to trade off performance versus complexity. As a rule of thumb, we 
have seen that Q — L/d max /A/J +1, where fd max is the maximum 
Doppler frequency and A/ is the subcarrier spacing, is an appropriate 
choice for Rayleigh fading [10]. In this paper, /d mox /A/ values are 
given as the normalized Doppler frequencies. It is concluded from 
these curves that the selection of the Q value is highly dependent on 
SNR values. In particular, for Eb/No = 20dB, different Q values 
show similar performance because the effects of ICI are not very 
obvious relative to the effects of the additive Gaussian noise. The 
Q value has a greater role for Eb/No above 25dB because ICI is 
dominant. We note that Q = 3 is sufficient for SNRs below 20dB. 

V. Conclusion 

Conventional detection methods such as ZF and MMSE have 
irreducible error floors at high normalized Doppler frequency fdT s 
since ICI corrupts the orthogonality among subcarriers. On the other 
hand, more sophisticated methods such as VBLAST require too much 
complexity, especially for large numbers of subcarriers. Therefore, 
we have proposed a new low-complexity maximum a posteriori 
probability (MAP) detection algorithm that provides excellent per- 
formance with manageable complexity for OFDM systems via the 
Gibbs sampling technique. In the simulation section, a comparison 
with other previously known receiver structures has been made and it 
has been demonstrated that MAP detection based on Gibbs sampling 



provides performance that is close to that of the optimal MAP 
detection algorithm for realistic fading conditions. 
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Abstract — In orthogonal frequency-division multiplexing (OFDM) 
systems operating over rapidly time-varying channels, the orthogonality 
between subcarriers is destroyed leading to inter-carrier interference 
(ICI) and resulting in an irreducible error floor. In this paper, a new 
and low-complexity maximum a posteriori probability (MAP) detection 
algorithm is proposed for OFDM systems operating over rapidly 
time-varying multipath channels. The detection algorithm exploits 
the banded structure of the frequency-domain channel matrix whose 
bandwidth is a parameter to be adjusted according to the speed of the 
mobile terminal. Based on this assumption, the received signal vector 
is decomposed into reduced dimensional sub-observations in such a 
way that all components of the observation vector contributing to the 
symbol to be detected are included in the decomposed observation 
model. The data symbols are then detected by the MAP algorithm 
by means of a Markov chain Monte Carlo (MCMC) technique in an 
optimal and computationally efficient way. Computational complexity 
investigation as well as simulation results indicate that this algorithm 
has significant performance and complexity advantages over existing 
suboptimal detection and equalization algorithms proposed earlier in 
the literature. 

Index Terms- OFDM, MAP detection, Monte carlo technique, Gibbs 
sampling. Intercarrier interference, fast time-varying channels. 



I. Introduction 

Orthogonal frequency-division multiplexing (OFDM) has been 
shown to be an effective method to overcome inter-symbol inter- 
ference (ISI) caused by frequency-selective fading with a simple 
transceiver structure, and is consequently used in several existing 
wireless local and metropolitan area standards such as the IEEE 
802.11 and IEEE 802.16 families. IEEE 802.11 wireless LAN 
(WLAN) technology has become very popular for providing data 
services to Internet users although its overall design and feature set 
are not well suited for outdoor broadband wireless access (BWA) 
applications [?]. Therefore, IEEE 802.16 has been developed as a new 
standard for BWA applications [?]. Recently, the much-anticipated 
Worldwide Interoperability for Microwave Access (WiMAX) tech- 
nology was introduced to promote the 802.16 standards while in- 
troducing features to enable mobile broadband services at vehicular 
speeds beyond 120 km/h. 

OFDM eliminates ISI and simply uses a one-tap equalizer to 
compensate multiplicative channel distortion in quasi-static channels. 
However, in fading channels with very high mobility, the time 
variation of the channel over an OFDM symbol period results 
in a loss of subchannel orthogonality which leads to inter-carrier 
interference (ICI). A considerable amount of research on OFDM 
receivers for quasi-static fading has been conducted, but a major 
hindrance to such receivers is the lack of mobility support [?]. Since 
mobility support is widely considered to be one of the key features 
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in wireless communication systems, and in this case ICI degrades 
the performance of OFDM systems, OFDM transmission over very 
rapidly time varying multipath fading channels has been considered 
recently in a number of recent works [?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, 
?]. 

The techniques proposed in these works range from linear equal- 
izers, based on the zero-forcing (ZF) or the minimum mean-squared 
error (MMSE) criterion [?, ?, ?, ?, ?, ?, ?, ?, ?, ?], to nonlinear 
equalizers based on decision-feedback or ICI cancelation [?, ?, 
?, ?, ?, ?]. Also near maximum-likelihood approaches have been 
proposed [?]. It has been shown that nonlinear equalizers based 
on ICI cancelation generally outperform linear approaches [?, ?, ?, 
?]. However, linear equalizers still preserve their importance mainly 
because they are less complex. 

In [?], performances of Matched Filter (MF), Least Squares (LS), 
Minimum Mean Square Error (MMSE) and MMSE with Successive 
Detection (SD) techniques with optimal ordering have been investi- 
gated. However, since the number of subcarriers is usually very large 
in high speed wide-band wireless standards, even the linear MMSE 
equalizer considered in [?] demands very high computational load. 

The specific structure of the Doppler-induced ICI in OFDM 
systems operating over highly mobile channels presents a distinctive 
feature of limited support of the Doppler spread that can be exploited 
by the receiver. References [?, ?, ?, ?, ?] exploit the banded character 
of the frequency-domain channel matrix to reach a complexity that 
is only linear in the number of subcarriers. In a certain sense, 
the assumption of a banded frequency-domain channel matrix is a 
natural extension of the time-invariant channel case, in which the 
frequency-domain channel matrix is diagonal and hence banded with 
the smallest possible bandwidth. 

In [?], using the banded structure of the channel matrix, a simple 
frequency domain equalizer has been proposed that can compensate 
for the loss of subchannel orthogonality due to ICI. However, the 
detection performance of the technique degrades substantially, since 
the data to be detected cannot fully use the contributing observation 
elements. The work presented in [?] combined [?] and [?] to derive 
a recursive decision feedback equalizer receiver for ICI suppression. 
The iterative MMSE serial linear equalizer (SLE) of [?], which takes 
the banded structure of the channel matrix into account, seems to be 
one of the most promising approaches to compensate for ICI. Iterative 
MMSE is then applied to estimate frequency-domain symbols. In [?], 
a block MMSE equalizer for OFDM systems operating over time- 
varying channels is presented. By exploiting the banded structure of 
the frequency-domain channel matrix, the complexity of the resulting 
algorithm turns out to be smaller than that of [?]. 

In this paper, a new computationally feasible, maximum a pos- 
teriori probability (MAP)-based data symbol detection algorithm is 
proposed for OFDM systems operating over highly mobile channels, 
as an alternative to the existing suboptimal equalization/detection 
techniques summarized in the above paragraphs. The proposed de- 
tection algorithm exploits the banded structure of the frequency- 
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domain channel matrix whose bandwidth is a parameter to be adjusted 
according to the speed of the mobile terminal. This assumption 
enables us to decompose the main received signal vector into finite 
numbers of reduced-dimensional, sub-received signal vectors from 
which the data symbols can be detected by the MAP algorithm in 
an optimal and computationally efficient way. The decomposition 
is achieved in such a way that all the components of the received 
vector that contribute to the symbol to be detected are included in the 
decomposed observation model. Data symbols in each sub-received 
signal model are then detected successively by a MAP detection 
algorithm. To implement MAP symbol detection in a computationally 
efficient way, we employ a Markov Chain Monte Carlo (MCMC) 
technique based on Gibbs sampling, which is a powerful statistical 
signal processing tool to estimate a posteriori probability (APP) 
values. 

The resulting detection algorithm is compared with previously 
proposed algorithms in terms of both bit error rate (BER) and 
complexity requirements. Computational complexity investigation as 
well as simulation results indicate that our algorithm has significant 
performance and complexity advantages over the existing suboptimal 
detection and equalization algorithms. 

II. System Model 

Let us consider an OFDM system with N subcarriers and available 
bandwidth B — 1/T„ where T 3 is the sampling period. A given 
sampling period is divided into N subchannels by equal frequency 
spacing A/ = B/N. At the transmitter, information symbols are 
mapped into possibly complex-valued transmitted symbols according 
to the modulation format employed. The symbols are processed 
by an N— length Inverse Fast Fourier Transform (IFFT) block that 
transforms the data symbol sequence into the time domain. The time- 
domain signal is extended by a guard interval containing G samples 
whose length is chosen to be longer than the expected delay spread 
to avoid ISI. The guard interval includes a cyclically extended part 
of the OFDM block to avoid ICI. Hence, the complete OFDM block 
duration is P — N + G samples. The resulting signal is converted to 
an analog signal by a digital-to-analog (D/A) converter. After shaping 
with a low-pass filter (e.g. a raised-cosine filter) with bandwidth B, 
it is transmitted through the transmit antenna with the overall symbol 
duration of PT S . 

Let h(m, I) represent the ith path (multipath component) of the 
time-varying channel impulse response at time instant t = mT s . The 
discrete-time received signal can then be expressed as follows: 



/(m) = >J h(m, l)d(m — I) + w(m), 



(1) 



where the transmitted signal d(m) at discrete sampling time m is 
given by 

1 JV-l 

d(m) = — = > s(k)e 



j2-Kmk/N 



(2) 



L is the total number of paths of the frequency selective fading 
channel, and w(m) is additive white Gaussian noise (AWGN) 
with zero mean and variance J3{|ii;(m)| 2 } = <rj. The sequence 
s(k),k = 0, 1, • • • , N — 1, in (2) represent either quadrature- 
amplitude modulation (QAM) or phase-sift-keying (PSK) modulated 
data symbols with E{\s(k)\ 2 } = 1. 

At the receiver, after passing through the analog-to-digital (A/D) 
converter and removing the cyclic prefix (CP), a fast Fourier trans- 
form (FFT) is used to transform the data back into the frequency 
domain. Lastly, the binary data is obtained after demodulation and 
channel decoding. 



The fading channel coefficients h(m, I) can be modeled as zero- 
mean complex Gaussian random variables. Based on the wide- 
sense stationary uncorrelated scattering (WSSUS) assumption, the 
fading channel coefficients in different paths are uncorrelated with 
each other. However, these coefficients are correlated within each 
individual path and have a lakes Doppler power spectral density 
having an autocorrelation function given by 

E{h(m,l)h*(n,l)} = a 2 hl J (27vf d T s (m-n)), (3) 

where o\ l denotes the power of the channel coefficients of the 
Zth path. f d is the Doppler frequency in Hertz so that the term 
fdT s represents the normalized Doppler frequency of the channel 
coefficients. Jo(.) is the zeroth order Bessel function of the first kind. 
By using (??) in (??), the received signal can be written as 

j iV-l L-l 

y{m) = -j= s(k)J2 h (™,l)e 3 * « +w(m), (4) 

fc=0 1=0 

which upon defining the time-varying channel transfer function 



H{k,m) = 2_^ l=0 h(m,l)e~ 



j2nlk/N 



(5) 



becomes 



y {m) = -L V s(k)H(k, m)e j27rmk/N + w(m). (6) 

The FFT output at the k th subcarrier, after excluding the guard 
interval, can be expressed as 



Y (k) = --L V y{m)e- ]2wmk/N 

v m— 

= s(k)G(k,k) + I(k) + W(k), 



(7) 



where I(k) is ICI caused by the time-varying nature of the channel 
given as 



I(k) = Y, s(i)G(k,i) 



(8) 



G(k, i) in (??) represents the average frequency domain time- 
varying channel response, defined as 



G(k,i) 4 (1/7Y) J2 H(i,m)e j2 * 



m(i-k)/N 



(9) 



Similarly, the term G(k, k) = i H ( k , m ) in ( ?? ) represent 

the portion of the average frequency domain channel response at the 
kth subcarrier and W(k) denotes discrete Fourier transform of the 
white Gaussian noise w(m): 



1 JV_1 
W(k) = -= > w(m)e 



j27vmk/N 



(10) 



Because of the term I(k) in (??), there is an irreducible error 
floor even in the training sequences since pilot symbols are also 
corrupted by ICI, arising from the fact that the time-varying channel 
destroys the orthogonality between subcarriers. Therefore, channel 
estimation should be performed either jointly with data or before the 
FFT block in order to compensate for the ICI. Note that if the channel 
is static or very slowly time-varying, that is H(k,m) ~ H(k), then 
it can be easily shown that G(k, k) = H(k) and I(k) — for 
k = 0, 1, • • ■ , N — 1, resulting in a received signal at the output of 
the FFT processor corresponding to the kth OFDM symbol given by 

Y(k) = H(k)s(k)+W(k). 
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From (6) and (7), the FFT output received signal can be expressed 
in vector form as 

Y = Gs + W (11) 

where Y = [Y(0), Y(l), Y(N- 1)] T , s = [ S (0), s(l), s(N- 
1)] T and W = [W (0), W (1), W(N - 1)] T . For k,i = 
0, 1, • • ■ , N - 1, the (k, i)th element of the matrix G = [G(k, i)] G 
qNxn re p resen ti n g the time- varying channel is given by (9). 

Under the assumption that the channel matrix G in (11) is per- 
fectly known at the receiver, the maximum likelihood (ML) detector 
performs an exhaustive search over the entire set of signal vectors 
whose components are selected from the signal constellation formed 
by the modulation scheme chosen. Especially in IEEE 802.16 based 
systems, the length N of each OFDM symbol is very large; it can 
take values as large as iV = 1024 or even iV = 2048 especially for 
high mobility applications. In this case an exhaustive search for the 
ML solution would be very complex since the search space has an 
extremely large number of constellation points, (jSI^, where \S\ is 
the cardinality of the signal constellation). On the other hand, all of 
the lower complexity linear detectors given in Table 1 are suboptimal 
since they do not take into account the correlation of the components 
of the transformed noise and yield noise enhancement. Recently, a 
nonlinear recursive detection technique using the decision-feedback 
principle, namely the MMSE-SD algorithm (VBLAST), has been 
proposed [?]. The performance of VBLAST depends critically on 
the order in which the data vector components are processed. To 
minimize error propagation effects and to improve the detection of 
unreliable components, more reliable data vector components should 
be detected first. Therefore, the algorithm depends on calculation 
of the post-detection signal-to-interference-plus-noise ratio (SINR) 
based upon MMSE detection as a measure of reliability, and so the 
calculation of SINR is compulsory at each iteration. Therefore this 
algorithm is computationally intensive as a number of pseudo inverse 
operations need to be performed. Moreover, its complexity grows 
exponentially with the total number of subcarriers. 

TABLE I 
Linear Detection Methods 



Method 


Solution 


Matched Filter (MF) 


s = Q{GTY} 


Zero Forcing (ZF) 


s = Q{(GTG)- 1 GTY} 


MMSE 


§ = Q{(G T G + I n <t* ) _1 G T Y} 



where Ijv is the N-by-N identity matrix. 



III. MAP Detection Algorithm 

The necessity of detecting large numbers of symbols in OFDM 
systems employed especially in highly mobile and wide-band wireless 
systems represents a significant computational burden as well as 
creating some convergence problems. However, as is known [?], 
time-varying channels produce a nearly-banded channel matrix G 
whose only main diagonal, Q subdiagonals and Q superdiagonals are 
nonzero. The bandwidth 2Q + 1, which is defined as the total number 
of non-zero diagonals in G, is a parameter to be adjusted according 
to the mobility-rate of the channel. The significant coefficients of 
G thus are confined to the 2Q + 1 central diagonals. The parameter 
Q G {0, 1, • • • N/2— 1} controls the target ICI response length; larger 
Q corresponds to a longer ICI span and, thus, increased estimation 
complexity. In general, Q should be chosen proportional to the width 
of the Doppler spectrum of the channel. 

We now present an optimal low-complexity MAP detection algo- 
rithm to detect the data symbols s from Y taking into account the 
banded structure of the channel matrix G. From the observation Y 



in (11), the receiver attempt to detect the OFDM output symbol s, 
assuming that G is completely known by the receiver. The banded 
structure of the channel matrix G implies that the data symbol 
s(k), k = 0, 1, • • • , N — 1, contributes to a maximum of 2Q + 1 
observation elements as follows: 

Y fc = [Y(j k ),Y(jk + !),••• ,Y{i h )f, for (12) 

jk = max{0, k — Q} and i k = min{N — 1, k + Q}. 

Based on this observation, the received signal in (11) can be de- 
composed into N reduced-dimensional sub-observations from which 
the data symbols can be detected in an optimal and computationally 
efficient way. For a given index k = 0, 1, • • • , N — 1 and Q, it can 
be easily shown from (11) and (12) that 

Y fc = G feSfc + W fc (13) 
where, s fe = [s(j fc )], W fe = [W(j k )], G fc = [G(i k ,j k )], for 

I Lk 4 max{0, k - Q} < i k < I Uk ± m in{7V -l,k + Q} (14) 
and 

J Lk = max{0,fc-2Q} < j k < J Uk = mva{N- 1, k+2Q}. (15) 

Note that due to the banded structure of G, some elements of the 
matrices G fe are zero and dim(G fc ) < (2Q + 1) x (2(2Q + 1) - 1). 
The Gfe's reach their maximum dimension when 2Q + 1 < k < 
(N-(2Q + 1)). 

For k — 0, 1, • • • , N — 1, the MAP estimate of the data symbol 
s(k) given Yfe is 

'sMAp(k) = s(k) = arg max P(s(k)\Y k ), (16) 

where S denotes the set of signal constellation points from which 
s(k) takes values. Based on this approach, s(k) can be detected 
sequentially for k = 0, 1, • • ■ , N — 1, incorporating the outcomes of 
the previous estimates in a decision-feedback mode as follows. 

■ For k = 0, determine the estimate S"(0) from (16). 

■ For k = k + 1 modify the observation vector Y^ by subtracting 
the terms coming from the contributions of the estimated data 
symbols S"(0),S"(1), • ■ • ,~s(k — 1) as 

fc-i 

Y fc 4 Yt - S^O ( 17 ) 

where is the Zth column of G^ and Jh k and Ju k are defined in 
(15). 

■ Determine the MAP estimate of ~s(k) from 

Y fc = G fc i fc + W fc (18) 

as 

s(k) = arg max P(s(k)\Y k ), (19) 

s(fc)es 

where 

6^G k -[gf> lg »V--,ri M ], (20 
and Sfc is the vector obtained by removing the first k — 1 elements 

Of Sfe. 

■ END IF k = N - 1. 

The major problem is finding the values of §MAp(k) in a com- 
putationally efficient manner. To see this difficulty we assume that 
the data symbols are independent and identically distributed binary 
phase shift keying (BPSK), taking values of +1 and — 1, Note that 
higher dimensional signal constellations can be treated similarly with 
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a straightforward extension. The conditional probability of s(k) given 
the observation vector Yk can be expressed as 

P(s(k) = +l\Y k ) =J2P(s(k) = +l,%|Y fe ) 

= P (*( k ) = +!|% Y fc )P(%| %) (21) 

where the second identity follows by applying the chain rule of 
probability. The vector in (21) is obtained by canceling the 
component s(k) in sfk and the summation is over all possible values 
of s^. The number of combinations that takes grows exponentially 
with the dimension of Sk and thus becomes prohibitive for large 
values of the size of this vector. Thus, we resort to the Gibbs sampler, 
a Monte Carlo method to calculate the a posteriori probabilities of 
the unknown symbols. 

A. MAP Detection Based on Gibbs Sampling 

The Gibbs sampler is an MCMC sampling method for numerical 
evaluation of multidimensional integrals. Its popularity is gained from 
the facts that it is capable of carrying out many complex Bayesian 
computations. In this section we briefly explain the application of the 
Gibbs sampling technique to our symbol detection problem where the 
observation process is given by (18). For notational convenience we 
drop the index k and the "tilde" from all the involved variables, e.g. 
Y, G are shorthand notations for Y^, G^, respectively. (21) can then 
be expressed as 

P(a(*) = +1|Y) = ]TP( S (&) = +l|s ¥ ,Y)P( Sr |Y) 

= £ S¥ | Y {P( S (fc) = +l|s ¥ ,Y)}. (22) 

According to the Gibbs sampling based statistical Monte Carlo 
estimation technique, an estimate of (22) can be evaluated by taking 
the empirical average 



P(s(k) = +1|Y) = ± ]T P(s(k) = +l|sjW Y) 



(23) 



where s~' for n = 1,2, ••• ,N e are samples drawn from the 
conditional distribution P(s^-|Y). There is a substantial body of 
literature concerning the Monte Carlo Gibbs sampling technique; see, 
e.g., [?], [?]. One possible version of the Gibbs sampler suitable for 
calculating the a posteriori probabilities in (21) may be summarized 
as follows. 

Let s = [s(0), s(l), • • • , s(N — 1)] T be a vector of unknown data 
symbols. Let Y be the observed signal. To generate random samples 
from the distribution P(s| Y), given the samples from the (n — l)th 



iteration s 



(n-1) 



„(«-!) 



(Oj.s^-^Cl),-.- ,s ( "- 1) (jV-l)]' r , the 
Gibbs algorithm iterates at the nth iteration as follows to generate 
the samples s (n) = [s<- n) (0), s (n) (l), ■ ■ ■ , s (n) (N - l)f: 

■ Initialize randomly; 

■ for n = 1, 2, • • • , N T and for k = 0, 1 • • • , N - 1, 

draw sample sj: n) from P ^s(fc)|sQ n) , • • • , 4-i> s L+i^> ' ' ' > s iv-i) 
Note that to ensure convergence, the Gibbs iteration is usually 
carried out for Nt = Nb + N s iterations. The first Nt iterations 
of the loop is called the bum-in period which is necessary for the 
Monte Carlo simulation to reach its stationary distribution. Only the 
samples s (n) = [s { n) , s[ n) , ■ ■ ■ , s^f, n = N b + 1, • • • , N T , from 
the last N s iterations are used to calculate the expectation. 
It is known that under regularity conditions [?, ?], 



(i) the distribution of s'™' converges to P(s|y), asn-> oo. 

(ii) (V^E"ii-P«fc) = +l|4 n) > Y ) = ^P(s(k) = 
+ 1|sf,Y)P(s ¥ |Y), asn^». 

B. Implementation of the Symbol Detector 

From the previous section, we recall that to compute P(s(k) = 
+l|Y fe ) = Y^s- P ( s ( k ) = + 1 > s fc|Yfc)> we need to perform the 
summation on the right-hand side of (22). When s has a large 
dimension, the exact evaluation of this summation may not be 
feasible and other more efficient techniques must be adopted. In this 
section the Gibbs sampling-based Monte Carlo method summarized 
in the previous section will be applied to develop a computationally 
efficient algorithm for calculation of the a posteriori probabilities 
P(s(fc)]Y). From (23), it follows that we need to evaluate P(s(k) = 
+l|s£ l) ,Y), for n = 1,2, ■■• ,N T . For this we define 



A 



(») 



P(s(k) = +l|s^,Y) 
= In ■ 



P(s{k) = 
from which it can be easily seen that 



-i|4"\y) 



P(s(k) = +l|s' n) ,Y 



i + «p(-aW) 



(24) 



(25) 



can be computed by expanding P(s(k) = +l|s^ l) , Y) as 



P(s(k) = +l|4 n) ,Y) = 



(26) 



p 


[Y\s(k) = +1, s<- n) )P( S (fc) = +1, s<- n) ) 




Es(t)6{+1 




(yk*) = +i,4*>) 


P( S (fc)=+l,S^) 



The data symbols are assumed to be independent and equally likely, 
Therefore, it follows from (26) and (24) that 



A 



(n) A 



hi 



p(YKfc) = +l,4 WJ ) 
p(Y|s(fe) = -l,sl n) ) 



(27) 



Since p(Y|s) ~ exp(— |y — Gs| 2 ), after some algebra, (27) can be 
expressed as 



A"( Sfc ) = -^5R{gi(Y-G ¥ s ¥ )} 

(Jill v J 



(28) 



where (-)t denotes the conjugate transpose and 5R{-} denotes the 
real part of its argument. is G with its fcth column gk removed. 
In summary, for k = 0, 1, • • • , N — 1, to estimate the a posteriori 
probabilities P(s(fe)|Y) in (21), the Gibbs sampler runs over all 

( In) ^ Nt 

symbols N a times to generate a collection of vectors \ s— > 

I k J n=JV b + l 

which are used in (23) to compute the desired quantities. 



C. Complexity Requirements 

The computational complexity of the MAP symbol detector based 
on Gibbs sampling proposed in this work is determined by the 
parameters N S ,Q,N and the constellation size of the transmitted data 
symbols. The computation of Yfe in (18) for k = 0, 1, • • ■ ,N — 1 
requires a maximum of (AQ 2 +2Q)N complex multiplications (CMs) 
and (4Q 2 +2Q)iV complex additions (CAs) per data block. Assuming 
BPSK signaling, the computation of the a posteriori probabilities 
in (22) requires a maximum of (8Q 2 + 2Q + l)NN s CMs and 
(8Q +2Q — 1)NN S CAs and computation of the empirical average 
of a posteriori probabilities of the data symbols in (23) requires NN S 
CSs. Therefore, the whole algorithm requires of maximum iV(4Q 2 + 
2Q + (4Q 2 + 2Q)N S ) CMs and A^(4Q 2 + 2Q + (4Q 2 + 2Q)N S ) 



